Joint Fairness Model with Applications to Risk Predictions for Under-represented Populations

In data collection for predictive modeling, under-representation of certain groups, based on gender, race/ethnicity, or age, may yield less-accurate predictions for these groups. Recently, this issue of fairness in predictions has attracted significant attention, as data-driven models are increasingly utilized to perform crucial decision-making tasks. Existing methods to achieve fairness in the machine learning literature typically build a single prediction model in a manner that encourages fair prediction performance for all groups. These approaches have two major limitations: i) fairness is often achieved by compromising accuracy for some groups; ii) the underlying relationship between dependent and independent variables may not be the same across groups. We propose a Joint Fairness Model (JFM) approach for logistic regression models for binary outcomes that estimates group-specific classifiers using a joint modeling objective function that incorporates fairness criteria for prediction. We introduce an Accelerated Smoothing Proximal Gradient Algorithm to solve the convex objective function, and present the key asymptotic properties of the JFM estimates. Through simulations, we demonstrate the efficacy of the JFM in achieving good prediction performance and across-group parity, in comparison with the single fairness model, group-separate model, and group-ignorant model, especially when the minority group's sample size is small. Finally, we demonstrate the utility of the JFM method in a real-world example to obtain fair risk predictions for under-represented older patients diagnosed with coronavirus disease 2019 (COVID-19).


Applied Context
The issue of making fair predictions has attracted significant attention recently in machine learning as a critical issue in the application of data-driven models. Though machine learning models are increasingly utilized to perform crucial decision-making tasks, recent evidence reveals that many carefully designed algorithms can learn biases from the underlying data and exploit these inequities when making predictions. For example, large systematic biases in prediction performance have been detected for machine learning models in areas such as recidivism prediction relative to race (Angwin et al., 2016), ranking of job candidates relative to gender (Lahoti et al., 2018) and face recognition relative to both race and gender (Ryu et al., 2017;Buolamwini and Gebru, 2018).
There is an emerging recognition that such biases in data often leads to unfair performance from predictive models in healthcare for certain groups (Char et al., 2018), such as women (Larrazabal et al., 2020), ethnic and racial minorities (Seyyed-Kalantari et al., 2021;Chen et al., 2021), and individuals with public insurance (Seyyed-Kalantari et al., 2021;Chen et al., 2021). Biased representations of different populations in biomedical studies, and the under-performance for under-represented populations, limit the potential benefits that can be achieved for these communities. In particular, when model-based predictions are used to prioritize patients for rationed services (e.g. organ transplantation, care management programs, or ICU services), under-performance for the under-represented patient populations will lead to unfair treatments for these patients (Paulus and Kent, 2020).
A particularly motivating example that we use in this paper is the mortality prediction for patients infected with coronavirus disease 2019 . As of January 23 2021, COVID-19 has infected more than 96 million people globally, accounting for more than 2 million known deaths. Older patients are particularly vulnerable to severe outcomes and death due to COVID-19. The Centers for Disease Control and Prevention (CDC) reported that the fatality rate was 18.8% for patients older than 80 years whereas the overall fatality rate is estimated at 5% or less for all patients (Kompaniyets and Goodman, 2021). This difference in survival highlights an urgent need for risk stratification of older patients with COVID-19 based on routine clinical assessments. However, most COVID-19 studies have not been stratified by age groups (Tehrani et al., 2021). Thus, as an example, when applying a risk prediction equation generated from the general population to older patients with COVID-19, the model in Tehrani et al. (2021) predicts high-risk scores overall due to their older age, higher prevalence of comorbidities and more laboratory abnormalities. This resulted in insufficient and unfair risk stratification for these patients as not all older patients are at the same risk of death from COVID-19 (Tehrani et al., 2021).

Existing Approaches
Methods to address fairness in the machine learning literature typically begin with a formal probabilistic definition of fairness. In the context of risk prediction, predictive fairness at the group level means that a risk prediction model has performance characteristics (such as accuracy, ranking, or calibration) that are relatively independent of group membership.
For example, if the false positive rate for a classification model is defined as P (ŷ = 1|y = 0), whereŷ is the model's prediction, then enforcing equality can be stated as requiring that these distributions be as close as possible between groups. Other definitions include demographic parity (Calders et al., 2009), equalized odds or equal opportunity (Hardt et al., 2016), disparate treatment, and impact and mistreatment (Zafar et al., 2019(Zafar et al., , 2017a. It is well-recognized that there is no unique optimal way to define fairness, leading to trade-offs between different approaches (Zafar et al., 2017b).
Given a fairness criterion, the second component of a fairness strategy requires an algorithmic approach, typically consisting of either 1) pre-processing the data by mapping the training data to a transformed space where the dependencies between sensitive attributes and class labels disappear (Kamiran and Calders, 2012;Dwork et al., 2018); or 2) postprocessing of a trained prediction model; for example, (Kamishima et al., 2012;Hardt et al., 2016) modify the probability of the decision being positive and negative predictions from an existing classifier to limit unfair discrimination; or 3) "in-process" methods, where fairness is accounted for during training of a model by adding a fairness constraint to the objective function. Examples of in-process methods include Zemel et al. (2013) who proposed to learn a fair representation of the data and classifier parameters by optimizing a non-convex function, and Zafar et al. (2017b) who defined a convex function as a measure of (un)fairness and suggested optimizing accuracy subject to the convex fairness constraints as well as their converse.
A key feature of nearly all existing approaches is that a single set of classifier parameters is estimated using fairness criteria that encourage fair prediction performance across all groups. This approach has two main limitations: i) fairness is often achieved by compromising accuracy of some groups; ii) the underlying relationship between dependent and independent variables may not be the same across groups, and the differences in predictive features may be of interest. In the example of predicting mortality risk for patients with COVID-19, while one would expect some features to have the same association with mortality for both older and younger patients, the associations between mortality and other features may be different between age groups. For instance, being overweight or obese (Body Mass Index [BMI] > 25kg/m 2 ) increases the risk for mortality associated with COVID-19, particularly among adults aged < 65 years (Kompaniyets and Goodman, 2021) However, geriatric BMI guidelines are different from those for younger adults. For older adults, higher BMIs are often associated with greater energy stores and a better nutritional state overall, which is beneficial for patients' survival outcomes when serious infections develop.
Estimating separate prediction models for each group does not leverage potential similarities between the groups. Moreover, estimating a single prediction model, even while using fairness criteria, will likely result in sub-optimal estimation or prediction performance for one group in order to achieve fair performance with a single set of parameters shared across groups. Outside the context of algorithmic fairness, Danaher et al. (2014) proposed the joint graphical lasso method, a technique for jointly estimating multiple models corresponding to distinct but related conditions. Their approach is based upon a penalized log-likelihood approach, which penalizes the differences between parameter estimates across groups. Penalized log-likelihood approaches have often been used by other authors like Lin (2007), Friedman et al. (2007b) etc. for similar estimation purposes while minimizing the disparities in estimates across groups. In all such cases, however, prediction performance was not emphasized.
In this paper, we propose the joint fairness model, a technique for jointly estimating multiple logistic regression models corresponding to distinct but related groups, in order to achieve fair prediction performance across groups. The model parameters are estimated by encouraging prediction fairness, while simultaneously ensuring high predictive accuracy irrespective of heterogeneity across groups. The rest of this paper is organized as follows. In Section 2, we present the proposed joint fairness model. Section 3 describes the algorithm to find its optimal solution. In Section 4, we discuss asymptotic consistency of the estimators.
We illustrate the performance of our proposed approach in simulation studies in Section 5; and apply our approach to the motivating example of predicting COVID-19 mortality outcomes for patients of different age groups in Section 6. Finally, we summarize and discuss our findings in Section 7.

Problem Formulation
For binary outcomes, consider K groups of datasets S k = {(X ki , y ki ) ∈ R p × {0, 1} : i = 1, · · · , n k } with K 2 representing group membership. Throughout the paper, group memberships are known and observed. Assuming that the n = K k=1 n k observations are independently distributed, then y ki ∼ Bernoulli(p ki ), andŷ ki : R p → {0, 1} is the predicted value based on predictor features X ki . We focus on the development of a fair prediction approach for the widely-used logistic regression model. The log-likelihood of the logistic model for the data from all groups takes the form Define β = (β 1 , · · · , β K ) ∈ R pK . Maximizing the likelihood function in (1) with respect to β k in each group separately yields the maximum likelihood estimates of parametersβ k for each group k, thus making separate predictionsŷ k per group. If we ignore group memberships,β can be estimated by maximizing the likelihood function in equation (1) setting all β k equal to a single global parameter vectorβ and making predictionsŷ per individual (irrespective of group) using that global parameter vector.
If the K datasets correspond to observations collected from K distinct but related groups, then one might wish to borrow strength across the K groups to estimate β and predictŷ, rather than estimating parameters β k for each group separately, or estimating a single set of β k for all k which could lead to heterogeneous prediction performance across the groups.
Therefore, instead of estimating β by maximizing the likelihood in equation (1), we consider a penalized log-likelihood approach and jointly estimate β by maximizing an objective function of K k=1 (β k ; X k , y k ) in equation (1) subject to constraints on (i) fairness, P F (β; X, y, λ F ) (ii) parameter similarity, P Sim (β; λ Sim ), and (iii) parameter sparsity, P Sp (β; λ Sp ).
We propose a fairness penalty function P F (β; X, y, λ F ) that encourages each group to have similar predictive performance. In this work we use equalized odds (Hardt et al., 2016) which encourages each group to have similar false positive rates (FPRs) and false negative rates (FNRs). Thus, we want to minimize the absolute difference between FPR j and FPR k |P (ŷ j = 1|y j = 0) − P (ŷ k = 1|y k = 0)|, and that between FNR j and FNR k : |P (ŷ j = 0|y j = 1) − P (ŷ k = 0|y k = 1)|. Under the logistic regression model, the absolute differences of FPRs and FNRs are nonconvex due to the nonconvexity of the sigmoid function.
We will instead minimize the absolute difference of the expected linear components of the two groups E X j β j y j = 0 − E X k β k y k = 0 . The proposition below shows that the absolute difference of the expected probabilities is upper-bounded by the absolute difference of the expected linear components. Thus, minimizing group differences in expected linear components guarantees that group difference of false predictions is minimized.
Note that the empirical estimate of the expectation is where S ky = {i : y ki = y} is a subgroup of subjects with a true response value y in group k, with y ∈ {0, 1}. Thus, our fairness penalty, that bridges the between-group gaps in the linear components of FPR k and FNR k , is defined as: where the summation j<k represents K k=2 k−1 j=1 for notational simplicity.
The similarity penalty P Sim (β; λ Sim ) is chosen to encourage similarity across the K parameter estimates. Here we use the generalized fused lasso penalty (Hoefling, 2010;Danaher et al., 2014;Dondelinger et al., 2018) herein referred to as the fusion-similarity penalty. The sparsity penalty P Sp (β) is chosen to encourage sparse estimates and to avoid ill-defined maximum likelihood estimates when sample sizes n k < p.
In the three penalty functions, λ F , λ Sim , and λ Sp are nonnegative hyperparameters. Here P Sp (β; λ Sp ), P F (β; X, y, λ F ), and P Sim (β; λ Sim ) are convex penalty functions, so that the objective function in equation (2) is convex in β. The proposed model jointly estimates β to achieve fair performance across groups, herein referred to as the joint fairness model (JFM).
In contrast, the dominant approach for fair predictions in the current literature is to estimate a single set of β parameters with constraints on the quality of performance metrics across groups. In the context of logistic regression, Bechavod and Ligett (2017) proposed a Single Fairness Model (SFM) to minimize the following objective function.
In contrast to the proposed JFM, it does not allow for group-specific different β k s, leading to two potential limitations: i) fairness is often achieved by compromising the likelihood ( (β; X k , y k )) of some groups; ii) inflexible model mis-specification when the β k s are different.
The proposed JFM objective function improves prediction parity through three components. First, it considers a weighted group total likelihood to upweight the groups with smaller sample sizes. Second, P F (β; X, y, λ F ) encourages estimates that achieve fair performance across groups. Third, the similarity penalty term improves estimation and prediction efficiencies when multiple subgroups are related but not identical (Hoefling, 2010;Danaher et al., 2014). Computationally, we found that multiple different combinations ofβ k often result in very similar objective values; thus, the similarity term can help optimization algorithms for the JFM converge to one of the multiple combinations by favoring similar values ofβ k . Other formats of the similarity penalty could be used in our proposed JFM framework. For example, the group lasso penalty (Yuan and Lin, 2006) has been shown to encourage similar sparsity patterns across groups (Obozinski et al., 2010;Danaher et al., 2014), while the fused lasso term is more aggressive in encouraging similarβ k estimates. The penalty functions in (3), (4), and (5) are based on the L 1 norm. They can be flexibly adapted to an L 2 penalization or to a combination of L 1 and L 2 penalizations. The differences between L 1 and L 2 penalties have been well discussed in the literature (Tibshirani, 1996;Zou and Hastie, 2005). For the fairness penalty, Bechavod and Ligett (2017) showed that there are no significant differences in the empirical predictive performances between the L 1 and L 2 fairness penalty forms other than the L 2 form of the similarity penalty penalizes large differences more aggressively so that models have less chance to obtain group-specific estimates.

Accelerated Smoothing Proximal Gradient Algorithm for JFM
In this section, we introduce an accelerated smoothing proximal gradient (ASPG) algorithm (Chen et al., 2012) to solve Problem (2) for the JFM. The objective function of (2) is convex in β so that a global optimal solution can be attained. However, conventional proximal gradientbased or coordinate descent approaches (generally used for lasso-like methods) cannot be directly applied to solve Problem (2) because there is no closed-form solution for the proximal operator associated with P FPR and P FNR .

Nesterov smooth approximation
To overcome the difficulty originating from the non-differentiability of the fairness and similarity penalties, we decouple the terms into a linear combination of the decision variables via the dual norm and then apply the Nesterov smoothing approximation (Nesterov, 2005).
We start with matrix representations of the fairness penalty terms P FPR (β; X, y, Similarly, we use the matrix representation of the similarity penalty P Sim (β; λ Sim ) = λ Sim Fβ 1 with F defined as below.
Here,X ky = 1 |S ky | i∈S ky X i is the average predictor vector for group k with true outcome y, and I p is the p-dimensional identity matrix. The matrix form of the fairness penalty term and the similarity penalty term is therefore defined as: Thus, the objective function (2) can be written in matrix form: where the associated proximal operator of D λ F ,λ Sim β 1 does not have a closed form solution.
We apply the Nesterov smooth approximation to approximate D λ F ,λ Sim β 1 by a smooth Proposition A.1 in Web Appendix 1 provides the maximum gap between D λ F ,λ Sim β 1 and its Nesterov approximation f µ (β; λ F , λ Sim ).
As demonstrated in Proposition A.2 in Web Appendix 1, for any µ > 0, f µ (β; λ F , λ Sim ) is smooth and convex with respect to β, whose gradient takes the following form: where where · 2 denotes the matrix spectral norm (which is equivalent to the largest singular value of the matrix). We can show further that α * can be calculated as where I is the indicator function. Details are provided in Web Appendix 1.

Accelerated Smoothing Proximal Gradient Algorithm
With whose first two terms are convex smooth functions. Although the sparsity penalty term k λ Sp k β k 1 is non-differentiable, it can be managed through the proximal gradient method using the soft-thresholding operator S with a closed form solution (Friedman et al., 2007a).
[Algorithm 1 about here.] Algorithm 1 presents the proposed ASPG algorithm, starting from parameter initialization, to gradient descent iterations with proximal and momentum steps, until convergence. Although Algorithm 1 minimizes the Nesterov smooth approximation, Theorem 3.1 proves that the solution can reach arbitrarily close to the global optimum of Problem (2).
Proof: See Web Appendix 4.
Based on Theorem 3.1, the rate of the convergence of Algorithm 1 is O pK δ(ε−δ) , given a desired accuracy ε > 0. The complexity of a single iteration of Algorithm 1 is O((n+K 2 )pK).
For additional details see Propositions A.4 and A.5 in Web Appendix 1. The proposed JFM Algorithm 1 is for the fusion-similarity term. The algorithm can be readily extended to include a group-similarity term as presented in Web Appendix 2.

Asymptotic properties of the JFM estimates
We now present asymptotic results for the JFM parameter estimatesβ obtained by solving Problem (2). We assume p remains constant and n increases to infinity.
Consider the following assumptions: Assumption 1. I(β k )/n k → C k , where C k is a positive definite p×p matrix, for k = 1, · · · , K, where I(β k ) is the information matrix of size p × p. For simplicity, we assume there are no intercept terms in β k .
Assumption 2. As n = min where I(β k ) is the empirical information matrix, and I p is a p × p identity matrix.
The following theorem proves √ n-consistency for the estimators, complying with the fairness and similarity constraints between the two groups as well as the sparsity constraint.
We note that the theorem holds even if the sample size of one group increases faster than the other group's.
Proof: See Web Appendix 4.

Simulation Study
We performed a series of simulations to evaluate the proposed JFM, and compared it with the approaches of a group-separate individual logistic regression model, a group-ignorant vanilla logistic regression model, and the SFM method (Bechavod and Ligett, 2017) implemented using an SFM-ASPG algorithm (see Web Appendix 3). When applying the group-separate model, regression coefficients were estimated for each group separately with an L 1 penalty.
The group-ignorant model estimates a single logistic regression with group membership as an additional covariate with an L 1 penalty.

Simulation Setup
We consider a two-group problem (K = 2) for simplicity, with group 1 as the over-represented group and group 2 as the under-represented group with respect to the sample sizes. The training samples were simulated as follows. The predictor matrix X k was independently generated from a standard normal distribution. The binary outcome y ki was then simulated from Bernoulli(π i (x ki )), where π i (x ki ) = exp(x ki β k ) 1+exp(x ki β k ) . Out of the total number of features, 40% in each group had non-zero coefficients (β's). The non-zero coefficients were each set to the value 3. The simulations were conducted under three scenarios to investigate performances at various levels of shared parameters, sample sizes and dimensionalities.
• Scenario 1 (Difference in True Model): The proportion of shared features between the two groups ranged from 0% to 100% of features with non-zero coefficients. The intercepts were selected so that the baseline event prevalences were at 30% and 50% for the underand over-represented groups The sample sizes were set at 500 and 200 for group 1 and 2 respectively. The number of features was set to p = 100. Note that when small proportions of features are shared between groups, the fairness and similarity penalty terms are misspecified given that Xβ k and β k are different between groups. As a result, performances in these settings allow us to investigate the robustness of the JFM approach to misspecification of the fairness and similarity penalty terms.
• Scenario 2 (Difference in Sample Size): The sample size of the under-represented group (group 2) ranged from 50 to 300 while the sample size of group 1 was fixed at 500.
The number of features was set to p = 100. Half of the features with non-zero coefficients were shared between the groups.
• Scenario 3 (High Dimensionality): The number of features p ranged from 50 to 2,000.
Sample sizes were 500 and 200 for group 1 and 2 respectively. For each value of p, 40 features had non-zero coefficients, with half of the non-zero features being shared.
We evaluated the methods on independent test datasets with large sample sizes (n = 1, 000 for both groups) under the same simulation setups. The Area under the Receiver Operating Characteristic curve (AUC) was used to assess the predictive ability of each model. Prediction unfairness was assessed by the group difference in AUCs. In addition, Mean Squared Errors (MSEs) for allβs were examined to assess parameter estimation performance, and TPRs of the associated features and TNRs of the non-associated features were used to assess variable selection performance. Medians and interquartile ranges (IQRs) of the assessment metrics were generated from 20 replicates for each experiment. Predictive performances and their unfairness in terms of FPR and FNR were calculated with a predicted probability cutoff of 0.5, as presented in Web Appendix 8. We present additional simulation scenarios in Web Appendix 9.

Choice of the Hyperparameters
The group-ignorant model, group-separate model, SFM, and JFM contain 1, K, 2, and K +2 hyperparameters respectively. For every method, five-fold cross-validation on the training dataset was used to determine the hyperparameters. For the vanilla models (group-separate and group-ignorant), the lasso penalty term was selected by optimizing cross-validation AUCs. For the fairness-aware models (SFM and JFM), we investigated and compared a series of evaluation metrics for selecting hyperparameters using cross-validation. The metrics include overall AUCs and Brier scores (defined as (p ki − y ki ) 2 ) on all samples (ignoring group memberships), group average of AUCs or Brier scores, and the group average of AUCs or Brier scores subtracting their disparities. Web Appendix 7, Web Figures 2, 3, and 4 demonstrate the empirical performance of the hyperparameters selected by the various strategies in the test datasets for the simulation scenarios. In summary, we find that the hyperparameters optimizing cross-validated group-average metrics showed better performance than sampleaverage or average metrics subtracting disparities. The hyperparameters optimizing AUCs in general generated the best AUC performances, while the hyperparameters optimizing Brier scores generated the best MSEs from the perspective of parameter estimation. We find that both are better empirically than threshold-based metrics such as classification accuracy. Lastly, the group average calculated by the harmonic mean is more robust than the arithmetic mean when the group sample sizes are unbalanced.

Simulation Results
For Scenario 1, Figure 1   The JFM showed consistently higher AUCs and smaller variances than those from all the other models. JFM outperforms the other models the most when the minority group's sample size is small, showing the benefits of borrowing information between groups in situations with unbalanced sample sizes. Figure 2(b) illustrates that the AUC of the majority group was not impacted for the JFM and group-separate methods. However, the AUC of the majority group decreased as the sample size of the under-represented group increased for the SFM and the group-ignorant models. This decrease highlights an undesirable aspect of these two methods, namely, they compromise accuracy by estimating a single set of classifier parameters. Figure   2(c) and 2(d) illustrates that the JFM achieves overall satisfactory AUCs and parity between groups across varying sample sizes of the under-represented group. Web Figure 6 compares the average of TPR and TNR and disparity of TPR and TNR of the four methods. In addition, the JFM substantially reduced coefficient MSEs (Figure 4(c)) under all simulated sample sizes for the minority group compared to all competing methods. The MSEs were similar between the JFM and the group separate model for the majority group, and both were lower than the MSEs of the group ignorant and SFM models. We investigated the empirical computational complexity of the JFM approach as a function of increasing numbers of features and sample sizes in Web Appendix 3. Web Figure 1 shows that the JFM computation time is approximately O(p 1.5 ) and O(n). This matches our theoretical analysis on complexity since the per-iteration complexity is O(np) and the rate of convergence is O(p 0.5 ). Details are presented in Web Appendix 6. We also implemented an alternative JFM algorithm using a group lasso similarity penalty (referred to as JFM-Group, Web Figure 15), and compared its performance with JFM-Fusion (presented above).
The results showed largely similar performances. The JFM-Fusion showed slightly better predictive performance than the JFM-Group when the two groups share more than 80% of the features, which is in line with previous reports that the fusion penalty term more aggressively enforces similarities across groups (Yuan and Lin, 2006;Obozinski et al., 2010;Danaher et al., 2014).

COVID-19 Risk Prediction Case Study
We applied the JFM, in comparison with other methods, to predict mortality related to COVID-19 from patients' routine ambulatory encounters and laboratory records prior to Fibrillation (AF); and routinely collected laboratory markers, such as lipid panels, blood panels, albumin, creatinine, aspartate aminotransferase (AST) etc. obtained from patients routine ambulatory histories before their COVID-19 infections. To build the prediction models, we randomly split the dataset into training (n = 8, 115, 70%) and testing (n = 3, 479, 30%) sets. We first standardized all features to zero-mean and unit variance. Five-fold crossvalidation was conducted on the training set to determine the hyperparameters for each model. Hyperparameters for the group-separate and group-ignorant models were selected to maximize the groupwise AUCs and the overall AUC, respectively, while those for the SFM and JFM were determined to maximize the harmonic mean of groupwise AUCs. Subsequently, we trained the final models with the optimal hyperparameters using the entire training set and applied the final models to the testing dataset to demonstrate their predictive performance. We repeated the training/testing split 10 times and averaged the performances across the 10 splits. Table 1 presents the AUCs and the averages of TPR and TNR of the four methods for each age group. The JFM performed better across all age groups than the separate model did, demonstrating that joint modeling yields higher efficiency. Compared with the group ignorant model, the JFM performed better in the three older age groups, with comparable AUC for the 50-64 age group, which resulted in smaller disparities in prediction performance overall. This phenomenon supports the observed pattern in simulation studies that the JFM reduced disparities in prediction performances without impacting those from the majority groups. In contrast, the SFM tended to reduce prediction disparities by lowering the performances for the majority groups. Figure 5 presents the boxplots of odds ratios (ORs) of selected demographic and clinical features estimated by the JFM. These results support the hypothesis that some features have common associations between groups, and some have group-specific ORs. For example, the decreasing OR estimates of BMI along agegroups confirmed the prior hypothesis that the association between BMI and COVID-19 mortality is heterogeneous between age-groups. In the JFM estimates, BMI is positively associated with higher risks of COVID-19 mortality for patients younger than 75, but with smaller and even reversed ORs in the oldest age groups. For older adults, higher BMIs are often associated with greater energy stores and a better nutritional state overall, which is beneficial for patients' survival outcomes when infected by COVID-19. The proportion of underweight patients (BMI<18) increased from 0.6% in the age group 50-64 to 5.5% in the age group 85+. The underweight status, often a proxy of frailness, has been repeatedly reported as a strong risk factor of COVID-19-induced multi-organ failure and mortality in older patients (Tehrani et al., 2021). On the other hand, the JFM can improve efficiencies for covariates with rare prevalence in a subgroup. For instance, dementia has been reported as a risk factor with COVID-19 mortality. In the group-separate model, dementia was insignificant in patients aged 50-64, mainly due to its low prevalence in this group (0.6%). In contrast, dementia was significantly associated with mortality in all age groups with similar ORs in the JFM estimates.

Conclusions and Discussion
In this study we introduced a joint fairness model for jointly estimating sparse parameters, on the basis of observations drawn from distinct but related groups, with the goal of achieving fair performances across groups. We employ an efficient accelerated smoothing proximal gradient algorithm to solve the joint fair objective function, which has convex penalty functions. Our algorithm is computationally tractable for high-dimensional datasets. Further, we presented the asymptotic distributions ofβ k . Our JFM predictions outperform competing approaches over a range of simulations and in an example application dataset.
We note that the JFM relies on separate hyperparameters (K + 2 hyperparameters) to control sparsity, fairness and similarity. This reliance can be viewed as a strength rather than a drawback because one can vary separately the amount of similarity, sparsity and fairness to enforce in the group-specific estimates. In situations with many groups, further assumptions can be made to reduce the number of sparsity hyperparameters (i.e. λ Sp k = c k λ Sp ). Possible choices of c k include 1 √ n k so that sparsity is inversely proportional to the number of samples.
While nearly all existing fairness-aware prediction approaches estimate a single set of classifier parameters, one exception is a recent study that proposes using multi-task learning (MTL) to improve algorithm fairness (Oneto et al., 2019). However, MTL research focuses on joint architecture, optimization, and task relationship learning, which is a different emphasis from the proposed JFM approach to improve risk prediction performance for underrepresented populations.
Theorem 4.1 established consistency of the JFM estimators. Its asymptotic distribution needs to be further investigated to lay the foundation of its inference. Moving forward, the proposed JFM framework can be extended for time-to-event outcomes by using similar constraints to those proposed here. It can also in principle be extended to non-linear models by adding a suitable fairness penalty term to the objective function.
Given the increasing ability to subclassify diseases according to their molecular features and the recognition that substantial heterogeneity exists in many molecular subtypes, most diseases will be eventually classified into a collection of multiple subtypes with unbalanced sample sizes. Therefore, the proposed JFM has wide application potential to improve prediction efficiencies and reduce subgroup prediction disparities beyond applications addressing gender, race/ethnicity and age disparities.

Data Availability Statement
The data that support the findings in this paper are available on request from the corresponding author. The data are not publicly available due to privacy or ethical restrictions.